Point Pattern Analysis: NASA Meteorite-Landings
Patrick Reza Schnurbusch
08/01/2021
The dataset contains the following variables:
name: the name of the meteorite (typically a location, often modified with a number, year, composition, etc) id: a unique identifier for the meteorite
nametype: one of:
- valid: a typical meteorite
- relict: a meteorite that has been highly degraded by weather on Earth
recclass: the class of the meteorite; one of a large number of classes based on physical, chemical, and other characteristics (see the Wikipedia article on meteorite classification for a primer)
mass: the mass of the meteorite, in grams
fall: whether the meteorite was seen falling, or was discovered after its impact; one of:
- Fell: the meteorite’s fall was observed
- Found: the meteorite’s fall was not observed
year: the year the meteorite fell, or the year it was found (depending on the value of fell)
reclat: the latitude of the meteorite’s landing
reclong: the longitude of the meteorite’s landing
GeoLocation: a parentheses-enclose, comma-separated tuple that combines reclat and reclong
Data Cleaning
We’ll clean out incorrect reccomended lat/long coordinates, we’ll also filter for USA via using a bounding box for lat/long within the most north/south/east/west bounds of the continental US.
meteorites <- meteorites %>%
dplyr::filter(year >= 860 & year <= 2016 ) %>% # filter out weird years
dplyr::filter(reclat != 0 | reclong != 0) %>%
dplyr::arrange(desc(year)) %>%
dplyr::filter(!is.na(mass)) %>%
dplyr::filter(mass > 0) %>%
dplyr::mutate(year = as.numeric(year))
# src for bbox
# https://gist.github.com/graydon/11198540
meteorites_clean <- meteorites %>%
dplyr::filter(reclat <= 49.3457868 & reclat >= 24.7433195) %>% #filter for north and south
dplyr::filter(reclong >= -124.7844079 & reclong <= -66.9513812) # filter for east and west Meteorites by Mass
pal <- colorNumeric(palette = "RdYlBu", domain = meteorites_clean$mass, reverse = T)
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng = ~reclong, lat = ~reclat,
color = ~pal(mass),
radius = 3,
stroke = FALSE, fillOpacity = 0.5) %>%
addLegend("bottomright", pal = pal, values = ~mass,
title = "Mass (grams)",
labFormat = labelFormat(suffix = " g"),
opacity = 1
)Heatmap of Meteorite impacts
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addHeatmap(lng = ~reclong, lat = ~reclat,
gradient = "viridis",
minOpacity = 0.1,
cellSize = 1,
radius = 3, blur = 15)Meteorites by Observed Status
Here we’ll visualize the meteorites via whether they were spotted on the ground after impact (“found”) or while in the mid-fall ("fell)
Overwhelmingly we can see most meteorites were spotted while on the ground after-impact was made, we have modified the size argument to help make the quantity more acute on our map.
pal <- colorFactor(c("#6EC5E9", "#FF1C51"), domain = c("Fell", "Found"))
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng = ~reclong, lat = ~reclat,
color = ~pal(fall),
radius = ~ifelse(fall == "Found", 1, 3),
stroke = FALSE, fillOpacity = 1) %>%
addLegend("bottomright", pal = pal, values = ~fall,
title = "Found Status",
opacity = 1)Distribution of Meteorites by Years
# Create a continuous palette function
pal <- colorBin(
palette = "Spectral",bins = 5,
domain = meteorites_clean$year, reverse = T)
leaflet(data = meteorites_clean, height=500, width=910) %>%
addFullscreenControl() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng = ~reclong, lat = ~reclat,
color = ~pal (year),
radius = 3,
stroke = FALSE, fillOpacity = 1) %>%
addLegend("bottomleft",
colors =c("#2B83BA", "#ABDDA4", "#FFFFBF","#FDAE61","#D7191C"),
labels = c("1600 - 1700", "1700 - 1800", "1800 - 1900", "1900 - 2000", "2000 - 2016"),
title = "Year Found",
opacity = 1)Most of the meteorites were discovered within last century (1900 - 2000)
Frequency distribtion by year
hchart(meteorites_clean$year,
color = "Purple",
name = "year") %>%
hc_yAxis(title = list(text = "Frequency (n)")) %>%
hc_xAxis(title = list(text = "<b> Year Observed </b>")) %>%
hc_add_theme(hc_theme_monokai())Density based analysis
Quadrat Density
We’ll divide the observation into a grid (h = 5, width = 10) and count the points that lie within each quadrat
# get the window for the data, i.e. the boundaries
W <- owin(xrange = c(min(meteorites_clean$reclong),
max(meteorites_clean$reclong)),
yrange = c(min(meteorites_clean$reclat),
max(meteorites_clean$reclat)))
# select just the lat/long columns
pp1 <- meteorites_clean %>%
dplyr::select(reclong,
reclat)
# convert to ppp object
meteorites_ppp <- as.ppp(pp1, W = W)
Q <- quadratcount(meteorites_ppp, nx= 5, ny=10)
# plot results
plot(meteorites_ppp, pch=20, cols="#FF1C51", main=NULL) # Plot points, with the tallies
plot(Q, add=TRUE) # Add quadrat grid
Computing the density
meteorites_ppp.km <- rescale(meteorites_ppp, 1000, "km")
# Compute the density for each quadrat (in counts per km2)
Q <- quadratcount(meteorites_ppp.km, nx= 6, ny=3)
Q.d <- intensity(Q)
# Plot the density
plot(intensity(Q, image=TRUE), main=NULL, las=1) # Plot density raster
plot(meteorites_ppp.km, pch=20, cex=0.6, col="#39FF14", add=TRUE) # Add points
Kernel density raster
K1 <- density(meteorites_ppp.km) # Using the default bandwidth
plot(K1, main=NULL, las=1)
contour(K1, add=TRUE)
Distance based analysis
We’ll apply three difference measures for distance bases analysis:
- The average nearest neighbor
- K & L function
- pair correlation function \(g\)
Average nearest neighbor analysis (ANN)
Computing the ANN distances between imapct sites. Compute the average first nearest neighbor distance (km), set k = 1
# set k - 1 for the first nearest neighbor distance
mean(nndist(meteorites_ppp.km, k = 1))## [1] 0.0002431465Compute the average second nearest neighbor, set k = 2
mean(nndist(meteorites_ppp.km, k = 2))## [1] 0.0003829143ANN Plot
ANN <- apply(nndist(meteorites_ppp.km, k=1:2000),2,FUN=mean)
# generate an ANN vs neighbor order plot
plot(ANN ~ eval(1:2000), type="b", main="ANN Plot", las=1)
The bottom axis shows the neighbor order number and the left axis shows the average distance (km). We can see a linear trend via
K & L Functions
Computing the K function
K <- Kest(meteorites_ppp.km)
plot(K, main=NULL, las=1)
This plot is visualizing the different estiamtes of \(K\) depending on the edge correction. The edge corrections plotted are listed in the legend.
\(K_{iso}\) refers to an isotropic implementation
\(K_{trans}\) refers to an translate implementation
\(K_{bord}\) refers to an border implementation
\(K_{pois}\) refers to an border implementation
Kernel density adjusted for covariate
Here we’ll estimate the meteorite point process intensity using the mass
In the following example, a Starbucks store point process’ intensity is estimated following the population density raster covariate
Hypothesis testing
Running an average nearest neighbor analysis for impact sites, in this we assume a uniform point density across space
ann.p <- mean(nndist(meteorites_ppp.km, k =1))
ann.p## [1] 0.0002431465The observed average nearest neighbor distance is 0.0002431465 km (0.24 m).
Generate the distribution of an expected
Hot Spot Spatial Analysis
Useful good here
Test for spatial autocorrelation/clustering in data
xy <- meteorites_clean[,c(9,8)]
meteorites_clean <- SpatialPointsDataFrame(coords = xy, data = meteorites_clean,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
print(nni(meteorites_clean))## Warning: data contain duplicated points## $NNI
## [1] 0.6192556
##
## $z.score
## [1] -30.13813
##
## $p
## [1] 1.534619e-199
##
## $expected.mean.distance
## [1] 0.3926433
##
## $observed.mean.distance
## [1] 0.2431465Our NNI score is less than 1, we can consider this as a good indication of clustering and can reject the Null Hypothesis that the points are uniform in pattern
Create the hotspot map
Define the hotspot map legend threshold
Areal Statistics
Meteorite Landings by State
In this section we’ll tally the number of meteorites via the State in which they land to gather some areal statistics.
# Downloading US State polygons
states <- states(cb = TRUE)##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|======== | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================= | 24%
|
|================= | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%# Convert to sf data
points <- st_as_sf(meteorites_clean)
# Set CRS
states <- st_transform(x = states, crs = 4326)
points <- st_transform(x = points, crs = 4326)
# Intersection between polygon and points ---------------------------------
intersection <- st_intersection(x = states, y = points)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries# Count number ofintersections per polygon (i.e. states)
states$pt_count <- lengths(st_intersects(states, points))
# Interactive Map
pal <- colorNumeric(
palette = "viridis",domain = states$pt_count, reverse = F, alpha = T)
labels_info = sprintf(
"<h2> %s </h2>
<hr>
Meteorites: %s ",
states$NAME,
states$pt_count
) %>% lapply(htmltools::HTML)
leaflet(states) %>%
addProviderTiles("CartoDB.DarkMatter") %>%
addPolygons(fillColor = ~pal(pt_count),
weight = 2,
opacity = 0.25,
color = "white",
dashArray = "1",
fillOpacity = 0.6,
highlight = highlightOptions(
weight = 3,
color = "#8ffcff",
fillColor = "#8ffcff",
dashArray = "1",
fillOpacity = 0.6,
bringToFront = TRUE),
label = labels_info,
labelOptions = labelOptions(
style = list("font-family" = "arial",
"font-style" = "bold",
"color" = "#FF1C51",
"background" = "black",
"border-color" = "#FFD300")
)
) %>%
setView(-98.5795, 39.8282, zoom=3) %>%
addLegend("bottomright",
pal = pal,
values = ~pt_count,
title = "Meteorites <hr>",
opacity = 1)